#############################################
# This file contains code to replicate      #
# Figures 1--3, A1, and A3--A7 in:          #
#                                           #
# "The Power to Hurt and the Effectiveness  #
# of International Sanctions"               #
#                                           #
# By Kerim Can Kavakli,                     #
# J. Tyson Chatagnier,                      #
# and Emre Hatipoglu                        #
#                                           #
# Forthcoming in the Journal of Politics    #
#                                           #
# File created 25 January 2019              #
#############################################

rm(list = ls())

library(readstata13)
library(ggplot2)
library(verification)
library(reshape2)

# setwd("")
                        # Uncomment this line to set a working directory

####################
# Source functions #
####################

source("auxiliary_functions.R")

############################
# Set seed for replication #
############################

set.seed(2018)

#############
# Load data #
#############

sanctions <- read.dta13("power2hurt_JOP_replication_dataset.dta")

sanctions <- sanctions[which(sanctions$trade_related == 1), ]
                        # Restrict the data set to only 
                        # trade-related sanctions

###############################
# Create a numeric Y variable #
###############################

sanctions$success <- ifelse(
                            sanctions$sanction_success == "successful sanction",
                            1,
                            0
                            )

######################
# Create an X matrix #
######################

x.mat <- cbind(
               sanctions$sender_mpower_over_targ,
               sanctions$targ_mpower_over_sender,
               sanctions$targ_gdppc,
               sanctions$targ_totalgdp,
               sanctions$targ_portfolio_concentration,
               sanctions$targ_cinc,
               sanctions$coalition_size,
               sanctions$targ_export_variety,
               sanctions$targ_portfolio_concentration * sanctions$targ_democracy,
               sanctions$targ_democracy,
               sanctions$coldwar,
               sanctions$security_issue,
               1
               )

colnames(x.mat) <- c(
                     "senderpower",
                     "targetpower",
                     "gdppc",
                     "gdptot",
                     "h_index",
                     "cap",
                     "num_sanctioners",
                     "num_ex",
                     "h_index_x_polity",
                     "polity",
                     "coldwar",
                     "security_issue",
                     "constant"
                     )


#################################
# Load results to plot Figure 1 #
#################################

beta <- read.delim(
                   "beta_model3.txt", 
                   header = TRUE
                   )

vcov <- read.delim(
                   "vcov_model3.txt", 
                   header    = TRUE, 
                   row.names = 1
                   )
                        # These files export the results from the logit
                        # analysis in Table 2, column 3 of the main text.

var.names <- c(
               "gdptot",
               "gdppc",
               "cap",
               "polity",
               "coldwar",
               "security_issue",
               "num_sanctioners",
               "senderpower",
               "targetpower",
               "num_ex",
               "h_index",
               "h_index_x_polity",
               "constant"
               )

beta <- as.numeric(beta)
vcov <- as.matrix(vcov)
                        # R reads these as data.frames, but we need
                        # them to be a vector and a matrix, respectively.

names(beta) <- var.names 
colnames(vcov) <- var.names
rownames(vcov) <- var.names 

###########################
# Reorder result matrices #
###########################

beta <- beta[colnames(x.mat)]

vcov <- vcov[colnames(x.mat), colnames(x.mat)]

##########################
# Set x profile as in LL #
##########################

profile.logit <- c(
                   apply(
                         x.mat[, 1 : 9],
                         2,
                         mean,
                         na.rm = TRUE
                         ),
                   apply(
                         x.mat[, (10) : 12],
                         2,
                         median,
                         na.rm = TRUE
                         ),
                   1
                  )
                        # Get means for continuous variables, set
                        # discrete variables to medians, and dichotomous
                        # variables to zero (incidentally, their medians
                        # in this case). We also include the constant.

profile.logit <- profile.logit[
                               order(
                                     match(
                                           names(profile.logit), 
                                           names(beta)
                                           )
                                     )
                               ]

############
# Figure 1 #
############

LogitPlots(
           x.mat        = x.mat,
           beta         = beta,
           vcov         = vcov,
           vals.profile = profile.logit,
           leg.pos      = c(
                            0.14,
                            0.175
                            )
           )

###############################
# Plot effects for threatened #
# and imposed sanctions from  #
# Modesl 4 and 5              #
###############################

beta.threat <- read.delim(
                          "beta_model4.txt", 
                          header = TRUE
                          )

vcov.threat <- read.delim(
                          "vcov_model4.txt", 
                          header    = TRUE, 
                          row.names = 1
                       )

beta.imp <- read.delim(
                       "beta_model5.txt", 
                       header = TRUE
                       )

vcov.imp <- read.delim(
                       "vcov_model5.txt", 
                       header    = TRUE, 
                       row.names = 1
                       )

var.names <- c(
               "gdptot",                      
               "gdppc",                      
               "cap",                      
               "polity",                      
               "coldwar",                      
               "security_issue",                      
               "num_sanctioners",                      
               "senderpower",                      
               "targetpower",                      
               "num_ex",                      
               "h_index",                      
               "h_index_x_polity",                      
               "constant"                      
               )

b.imp <- as.numeric(beta.imp)
vcov.imp <- as.matrix(vcov.imp)
b.threat <- as.numeric(beta.threat)
vcov.threat <- as.matrix(vcov.threat)
                        # R reads these as data.frames, but we need
                        # them to be a vector and a matrix, respectively.

names(b.imp) <- var.names 
names(b.threat) <- var.names 
colnames(vcov.imp) <- var.names
rownames(vcov.imp) <- var.names 
colnames(vcov.threat) <- var.names
rownames(vcov.threat) <- var.names 

###########################
# Reorder result matrices #
###########################

b.imp <- b.imp[colnames(x.mat)]
b.threat <- b.threat[colnames(x.mat)]
vcov.imp <- vcov.imp[colnames(x.mat), colnames(x.mat)]
vcov.threat <- vcov.threat[colnames(x.mat), colnames(x.mat)]

############
# Figure 2 #
############

LogitPlots(
           x.mat = x.mat,
           beta  = b.threat,
           vcov  = vcov.threat,
           vals.profile = profile.logit,
           leg.pos = c(
                       0.14,
                       0.825
                       )
           )

############
# Figure 3 #
############

LogitPlots(
           x.mat = x.mat,
           beta  = b.imp,
           vcov  = vcov.imp,
           vals.profile = profile.logit,
           leg.pos = c(
                       0.14,
                       0.825
                       )
           )

#############################
# Prepare data for Figure 4 #
#############################

x.ll <- cbind(
              sanctions$sender_mpower_over_targ,
              sanctions$targ_mpower_over_sender,
              sanctions$targ_gdppc,
              sanctions$targ_totalgdp,
              sanctions$targ_portfolio_concentration,
              sanctions$targ_cinc,
              sanctions$coalition_size,
              sanctions$targ_export_variety,
              sanctions$targ_democracy,
              sanctions$coldwar,
              sanctions$security_issue
              )

colnames(x.ll) <- c(
                    "senderpower",
                    "targetpower",
                    "gdppc",
                    "gdptot",
                    "h_index",
                    "cap",
                    "num_sanctioners",
                    "num_ex",
                    "polity",
                    "coldwar",
                    "security_issue"
                    )

######################
# Set q1, q2, and q3 #
######################

q1 <- 8
                        # Treating number of sanctioners and
                        # number of exporters as continuous
q2 <- 1
q3 <- 2

########################################
# Standardize continuous variables for #
# numerical stability.                 #
########################################

x.ll[, 1 : q1] <- apply(
                        x.ll[, 1 : q1],
                        2,
                        StandVar
                        )

##############################
# Get locations of variables #
##############################

loc.send <- which(colnames(x.ll) == "senderpower")
loc.targ <- which(colnames(x.ll) == "targetpower")
loc.num <- which(colnames(x.ll) == "num_ex")
loc.h <- which(colnames(x.ll) == "h_index")

########################
# Set profile for plot #
########################

vals.profile <- c(
                  apply(
                        x.ll[, 1 : q1],
                        2,
                        mean,
                        na.rm = TRUE
                        ),
                  apply(
                        x.ll[, (q1 + 1) : ncol(x.ll)],
                        2,
                        median,
                        na.rm = TRUE
                        )
                  )
                        # We will want to find average effects, setting 
                        # continuous variables to means and discrete 
                        # variables to medians

####################
# Cross-validation #
####################

# cv.grid <- expand.grid(
#                        c(
#                          4, 
#                          10, 
#                          25, 
#                          50
#                          ), 
#                        c(
#                          0.1, 
#                          0.25, 
#                          0.5, 
#                          0.75
#                          ), 
#                        c(
#                          0.1, 
#                          0.25, 
#                          0.5, 
#                          0.75
#                          )
#                        )
                        # Use this to run cross-validation on a grid of 64
                        # different starting values.

# cv <- crossValidation(
#                       stval.mat = cv.grid, 
#                       X         = x.ll, 
#                       Y         = sanctions$success, 
#                       q1        = q1, 
#                       q2        = q2, 
#                       q3        = q3,
#                       nclus     = 4
#                       )

# bw <- as.vector(cv$par[cv$llik == max(cv$llik), ])
                        # Uncomment the lines above to run the
                        # cross-validation from scratch. The
                        # next line should then be commented out.

bw <- c(
        7.0056673, 
        0.7984081, 
        0.5170194
        )
                        # The chosen smoothing parameters should be those
                        # that maximize the log-likelihood.
                        # The values selected by this procedure were:
                        # (7.0056673, 0.7984081, 0.5170194)
                        # With log-likelihood:
                        # -383.4029


############
# Figure 4 #
############

MakePlot(
         sanctions$success,
         bw,
         x.ll,
         leg.pos = c(
                     0.14,
                     0.175
                     )
         )

######################################
#                                    #
#          APPENDIX FIGURES          #
#                                    #
######################################

############################
# Plot substantive effects #
# for Figure A1            #
############################

x.a4 <- cbind(
              sanctions$sender_mpower_over_targ,
              sanctions$targ_mpower_over_sender,
              sanctions$targ_gdppc,
              sanctions$targ_totalgdp,
              sanctions$targ_portfolio_concentration,
              sanctions$targ_cinc,
              sanctions$coalition_size,
              sanctions$targ_export_variety,
              sanctions$targ_portfolio_concentration * sanctions$targ_democracy,
              sanctions$targ_democracy,
              sanctions$coldwar,
              sanctions$security_issue,
              sanctions$smart_sanction,
              sanctions$us_primarysender,
              sanctions$igo_involved,
              1
              )

colnames(x.a4) <- c(
                    "senderpower",
                    "targetpower",
                    "gdppc",
                    "gdptot",
                    "h_index",
                    "cap",
                    "num_sanctioners",
                    "num_ex",
                    "h_index_x_polity",
                    "polity",
                    "coldwar",
                    "security_issue",
                    "smart_sanction",
                    "us_primarysender",
                    "igo_involved",
                    "constant"
                    )

b.a4 <- read.delim(
                   "beta_a4_m1.txt", 
                   header = TRUE
                   )

vcov.a4 <- read.delim(
                      "vcov_a4_m1.txt", 
                      header    = TRUE, 
                      row.names = 1
                      )

var.names <- c(
               "gdptot",
               "gdppc",
               "cap",
               "polity",
               "coldwar",
               "security_issue",
               "num_sanctioners",
               "senderpower",
               "targetpower",
               "num_ex",
               "h_index",
               "h_index_x_polity",
               "smart_sanction",
               "us_primarysender",
               "igo_involved",
               "constant"
               )

b.a4 <- as.numeric(b.a4)
vcov.a4 <- as.matrix(vcov.a4)

names(b.a4) <- var.names 
colnames(vcov.a4) <- var.names
rownames(vcov.a4) <- var.names 

b.a4 <- b.a4[colnames(x.a4)]
vcov.a4 <- vcov.a4[colnames(x.a4), colnames(x.a4)]

############################
# Create profile of values #
############################

profile.a4 <- c(
                apply(
                      x.a4[, 1 : 9],
                      2,
                      mean,
                      na.rm = TRUE
                      ),
                apply(
                      x.a4[, 10 : 15],
                      2,
                      median,
                      na.rm = TRUE
                      ),
                1
                )
                        # Get means for continuous variables, set
                        # discrete variables to medians, and dichotomous
                        # variables to zero (incidentally, their medians
                        # in this case). We also include the constant.

profile.a4 <- profile.a4[
                         order(
                               match(
                                     names(profile.a4), 
                                     names(beta)
                                     )
                               )
                         ]


################
# Create plots #
################

orig.plots <- LogitPlots(
                         x.a4, 
                         b.a4, 
                         vcov.a4, 
                         profile.a4,
                         leg.pos = c(
                                     0.14,
                                     0.175
                                     ),
                         .plot = FALSE
                         )

p.dem <- GetPreds(
                  xmat        = x.a4,
                  coefs       = b.a4,
                  covmat      = vcov.a4,
                  change.lo   = -10,
                  change.hi   = 10,
                  var.change  = 10,
                  interaction = c(
                                  5, 
                                  9
                                  ),
                  app.fn      = profile.a4,
                  tval        = 1.645,
                  n.out       = 21,
                  ci          = TRUE
                  )

vals.h.lo <- profile.a4
vals.h.lo[5] <- min(
                    sanctions$targ_portfolio_concentration,
                    na.rm = TRUE
                    )

vals.h.hi <- profile.a4
vals.h.hi[5] <- max(
                    sanctions$targ_portfolio_concentration,
                    na.rm = TRUE
                    )

p.dem.lo <- GetPreds(
                     xmat        = x.a4,
                     coefs       = b.a4,
                     covmat      = vcov.a4,
                     change.lo   = -10,
                     change.hi   = 10,
                     var.change  = 10,
                     interaction = c(
                                     5, 
                                     9
                                     ),
                     app.fn      = vals.h.lo,
                     tval        = 1.645,
                     n.out       = 21,
                     ci          = TRUE
                     )

p.dem.hi <- GetPreds(
                     xmat        = x.a4,
                     coefs       = b.a4,
                     covmat      = vcov.a4,
                     change.lo   = -10,
                     change.hi   = 10,
                     var.change  = 10,
                     interaction = c(
                                     5, 
                                     9
                                     ),
                     app.fn      = vals.h.hi,
                     tval        = 1.645,
                     n.out       = 21,
                     ci          = TRUE
                     )

df.dem <- data.frame(
                     Change      = p.dem$change,
                     Hi          = p.dem$hi,
                     Lo          = p.dem$lo,
                     Probability = p.dem$pred
                     )

df.dem.lo <- data.frame(
                        Change      = p.dem.lo$change,
                        Hi          = p.dem.lo$hi,
                        Lo          = p.dem.lo$lo,
                        Probability = p.dem.lo$pred
                        )

df.dem.hi <- data.frame(
                        Change      = p.dem.hi$change,
                        Hi          = p.dem.hi$hi,
                        Lo          = p.dem.hi$lo,
                        Probability = p.dem.hi$pred
                        )
plot.dem <- ggplot(
                   df.dem,
                   aes(
                       x = Change,
                       y = Probability
                       ) 
                   ) +
                   geom_point(
                              size   = 1.5,
                              colour = "navyblue"
                              ) +
                   geom_errorbar(
                                 aes(ymin = Lo,
                                     ymax = Hi
                                     ),
                                 colour = "navyblue",
                                 size   = 0.75
                                 ) +
                   xlab("Polity Score") +
                   ylab("Probability of Sanctions Success") +
                   theme_bw()

plot.dem.lo <- ggplot(
                      df.dem.lo,
                      aes(
                          x = Change,
                          y = Probability
                          ) 
                      ) +
                      geom_point(
                                 size = 1.5,
                                 colour = "navyblue"
                                 ) +
                      geom_errorbar(
                                    aes(
                                        ymin = Lo,
                                        ymax = Hi
                                        ),
                                    colour = "navyblue",
                                    size   = 0.75
                                    ) +
                      xlab("Polity Score") +
                      ylab("Probability of Sanctions Success") +
                      theme_bw()

plot.dem.hi <- ggplot(
                      df.dem.hi,
                      aes(
                          x = Change,
                          y = Probability
                          ) 
                      ) +
                      geom_point(
                                 size   = 1.5,
                                 colour = "navyblue"
                                 ) +
                      geom_errorbar(
                                    aes(
                                        ymin = Lo,
                                        ymax = Hi
                                        ),
                                    colour = "navyblue",
                                    size   = 0.75
                                    ) +
                      xlab("Polity Score") +
                      ylab("Probability of Sanctions Success") +
                      theme_bw()

df.three <-data.frame(
                      Change      = rep(
                                        p.dem$change,
                                        times = 3
                                        ),
                      Hi          = c(
                                      p.dem.lo$hi,
                                      p.dem$hi,
                                      p.dem.hi$hi
                                      ),
                      Lo          = c(
                                      p.dem.lo$lo,
                                      p.dem$lo,
                                      p.dem.hi$lo
                                      ),
                      Probability = c(
                                      p.dem.lo$pred,
                                      p.dem$pred,
                                      p.dem.hi$pred
                                      ),
                      HH          = rep(
                                        c(
                                          "Low",
                                          "Mean",
                                          "High"
                                          ),
                                        each = length(p.dem$change)
                                        )
                      )

df.three$HH <- factor(
                      df.three$HH, 
                      c(
                        "Low", 
                        "Mean", 
                        "High"
                        )
                      )
                        # Reorder the values for plotting.

plot.dem.all <- ggplot(
                       df.three,
                       aes(
                           x        = Change,
                           y        = Probability,
                           colour   = HH,
                           fill     = HH,
                           linetype = HH
                           ) 
                       ) +
                       geom_line(size   = 1.10) +
                       geom_ribbon(
                                   aes(
                                       ymin = Lo,
                                       ymax = Hi,
                                       ),
                                   alpha  = 0.3,
                                   colour = NA
                                   ) +
                       xlab("Polity Score") +
                       ylab("Probability of Sanctions Success") +
                       scale_linetype_discrete(name = "Portfolio Concentration") +
                       scale_colour_discrete(name = "Portfolio Concentration") +
                       scale_fill_discrete(name = "Portfolio Concentration") +
                       ylim(
                            0, 
                            1
                            ) +
                       theme_bw() +
                       theme(
                             legend.position = c(
                                                 0.15,
                                                 0.85
                                                 )
                             ) 

########################
# And plot IGO effects #
########################


p.igo <- GetPreds(
                  xmat       = x.a4,
                  coefs      = b.a4,
                  covmat     = vcov.a4,
                  var.change = 15,
                  change.lo  = 0,
                  change.hi  = 1,
                  app.fn     = profile.a4,
                  n.out      = 2,
                  tval       = 1.645,
                  ci         = TRUE
                  )

df.igo <- data.frame(
                     Change      = c(
                                     "Not Inolved",
                                     "Involved"
                                     ),
                     Hi          = p.igo$hi,
                     Lo          = p.igo$lo,
                     Probability = p.igo$pred
                     )

plot.igo <- ggplot(
                   df.igo,
                   aes(
                       x = Change,
                       y = Probability
                       ) 
                   ) +
                   geom_point(
                              size   = 3.5,
                              colour = "navyblue"
                              ) +
                   geom_errorbar(
                                 aes(
                                     x    = Change,
                                     ymin = Lo,
                                     ymax = Hi
                                     ),
                                 colour = "navyblue",
                                 size   = 1.2,
                                 width  = 0.1
                                 ) +
                   xlab("IGO Involvement") +
                   ylab("Probability of Sanctions Success") +
                   ylim(
                        0, 
                        1
                        ) +
                   theme_bw()

#############
# Figure A1 #
#############

multiplot(
          orig.plots$send,
          orig.plots$h,
          plot.dem.all,
          orig.plots$targ,
          orig.plots$num,
          plot.igo,
          cols = 2
          )

###########################
# Read in status quo data #
# for Figure A2           #
###########################

data.sq <- read.dta13("sq_observations.dta")

##########################
# Combine with main data #
##########################

cols.diff <- c(
               1,
               3,
               9
               )

colnames(data.sq)[cols.diff] <- c(
                                  "targetstate",
                                  "primarysender",
                                  "target_polity2"
                                  )

sanctions.sq <- merge(
                      sanctions,
                      data.sq,
                      all = TRUE
                      )

x.sq <- cbind(
              sanctions.sq$sender_mpower_over_targ,
              sanctions.sq$targ_mpower_over_sender,
              sanctions.sq$targ_gdppc,
              sanctions.sq$targ_totalgdp,
              sanctions.sq$targ_portfolio_concentration,
              sanctions.sq$targ_cinc,
              sanctions.sq$targ_export_variety,
              sanctions.sq$targ_democracy,
              sanctions.sq$coldwar
              )

colnames(x.sq) <- c(
                    "senderpower",
                    "targetpower",
                    "gdppc",
                    "gdptot",
                    "h_index",
                    "cap",
                    "num_ex",
                    "polity",
                    "coldwar"
                    )
q1.sq <- 7
q2.sq <- 1
q3.sq <- 1

x.sq[, 1 : q1.sq] <- apply(
                        x.sq[, 1 : q1.sq],
                        2,
                        StandVar
                        )

loc.send <- which(colnames(x.sq) == "senderpower")
loc.targ <- which(colnames(x.sq) == "targetpower")
loc.num <- which(colnames(x.sq) == "num_ex")
loc.h <- which(colnames(x.sq) == "h_index")

vals.profile <- c(
                  apply(
                        x.sq[, 1 : q1.sq],
                        2,
                        mean,
                        na.rm = TRUE
                        ),
                  apply(
                        x.sq[, (q1.sq + 1) : ncol(x.sq)],
                        2,
                        median,
                        na.rm = TRUE
                        )
                  )

####################
# Cross-validation #
####################

# cv <- crossValidation(
#                       stval.mat = cv.grid, 
#                       X         = x.sq, 
#                       Y         = sanctions.sq$success, 
#                       q1        = q1.sq, 
#                       q2        = q2.sq, 
#                       q3        = q3.sq,
#                       nclus     = 4
#                       )

# bw <- as.vector(cv$par[cv$llik == max(cv$llik), ])
#                         # Uncomment the lines above to run the
#                         # cross-validation from scratch. The
#                         # next line should then be commented out.

bw.sq <- c(
           5.7827382, 
           0.8218152, 
           0.4192080
           )
                        # The chosen smoothing parameters should be those
                        # that maximize the log-likelihood.
                        # The values selected by this procedure were:
                        # (5.7827382, 0.8218152, 0.4192080)
                        # With log-likelihood:
                        # -378.8643

#############
# Figure A2 #
#############

MakePlot(
         sanctions.sq$success,
         bw.sq,
         x.sq,
         leg.pos = c(
                     0.14,
                     0.175
                     ),
         q1.val = q1.sq,
         q2.val = q2.sq,
         q3.val = q3.sq
         )

############################
# Create Figures A3 and A4 #
############################

vals.profile <- c(
                  apply(
                        x.ll[, 1 : q1],
                        2,
                        mean,
                        na.rm = TRUE
                        ),
                  apply(
                        x.ll[, (q1 + 1) : ncol(x.ll)],
                        2,
                        median,
                        na.rm = TRUE
                        )
                  )
                        # Recreate vals.profile with full data

sanctions$vic_threat <- ifelse(
                               sanctions$success == 1 &
                               sanctions$imposed == 0,
                               1,
                               0
                               )
                        # Sanctions that fail and sanction threats
                        # that lead to imposition are all considered
                        # failures here.

loc.send <- which(colnames(x.ll) == "senderpower")
loc.targ <- which(colnames(x.ll) == "targetpower")
loc.num <- which(colnames(x.ll) == "num_ex")
loc.h <- which(colnames(x.ll) == "h_index")

################################
# Cross-validation for threats #
################################

# cv <- crossValidation(
#                       stval.mat = cv.grid, 
#                       X         = x.mat, 
#                       Y         = sanctions$vic_threat, 
#                       q1        = q1, 
#                       q2        = q2, 
#                       q3        = q3,
#                       nclus     = 4
#                       )

# bw.threat <- as.vector(cv$par[cv$llik == max(cv$llik), ])
                        # Uncomment the lines above to run the
                        # cross-validation from scratch. The
                        # next line should then be commented out.

bw.threat <- c(
               10.7459074, 
               0.8119127, 
               0.4652445
               )
                        # The chosen smoothing parameters should be those
                        # that maximize the log-likelihood.
                        # The values selected by this procedure were:
                        # (10.7459074, 0.8119127, 0.4652445)
                        # With log-likelihood:
                        # -358.4863

#############
# Figure A3 #
#############

MakePlot(
         sanctions$vic_threat,
         bw.threat,
         x.ll,
         leg.pos = c(
                     0.14,
                     0.75
                     )
         )

##########################################
# Cross-validation for imposed sanctions #
##########################################

sanc.imp <- sanctions[sanctions$imposed == 1, ]

x.imp <- cbind(
               sanc.imp$sender_mpower_over_targ,
               sanc.imp$targ_mpower_over_sender,
               sanc.imp$targ_gdppc,
               sanc.imp$targ_totalgdp,
               sanc.imp$targ_portfolio_concentration,
               sanc.imp$targ_cinc,
               sanc.imp$coalition_size,
               sanc.imp$targ_export_variety,
               sanc.imp$targ_democracy,
               sanc.imp$coldwar,
               sanc.imp$security_issue
               )

colnames(x.imp) <- c(
                     "senderpower",
                     "targetpower",
                     "gdppc",
                     "gdptot",
                     "h_index",
                     "cap",
                     "num_sanctioners",
                     "num_ex",
                     "polity",
                     "coldwar",
                     "security_issue"
                     )

x.imp[, 1 : q1] <- apply(
                         x.imp[, 1 : q1],
                         2,
                         StandVar
                         )

vals.profile <- c(
                  apply(
                        x.imp[, 1 : q1],
                        2,
                        mean,
                        na.rm = TRUE
                        ),
                  apply(
                        x.imp[, (q1 + 1) : ncol(x.imp)],
                        2,
                        median,
                        na.rm = TRUE
                        )
                  )


# cv <- crossValidation(
#                       stval.mat = cv.grid, 
#                       X         = x.imp, 
#                       Y         = sanc.imp$success, 
#                       q1        = q1, 
#                       q2        = q2, 
#                       q3        = q3,
#                       nclus     = 4
#                       )

# bw <- as.vector(cv$par[cv$llik == max(cv$llik), ])
                        # Uncomment the lines above to run the
                        # cross-validation from scratch. The
                        # next line should then be commented out.

bw.imp <- c(
            15.7940815, 
            0.9006885, 
            0.2626998
            )
                        # The chosen smoothing parameters should be those
                        # that maximize the log-likelihood.
                        # The values selected by this procedure were:
                        # (15.7940815, 0.9006885, 0.2626998)
                        # With log-likelihood:
                        # -188.4712

#############
# Figure A4 #
#############

MakePlot(
         sanc.imp$success,
         bw.imp,
         x.imp,
         leg.pos = c(
                     0.14,
                     0.175
                     )
         )

######################
# Compare models for #
# Figures A5 and A6  #
######################

x.full <- cbind(
                sanctions$targ_totalgdp,
                sanctions$targ_gdppc,
                sanctions$targ_cinc,
                sanctions$targ_democracy,
                sanctions$coldwar,
                sanctions$security_issue,
                sanctions$coalition_size,
                sanctions$sender_mpower_over_targ,
                sanctions$targ_mpower_over_sender,
                sanctions$trade_dependence,
                sanctions$targ_export_variety,
                sanctions$targ_portfolio_concentration,
                sanctions$targ_portfolio_concentration * sanctions$targ_democracy,
                1
                )

x.pow <- x.full[, -10]
x.dep <- x.full[, -(8 : 9)]

b.pow <- as.numeric(
                    read.delim(
                               "beta_a9_m1.txt", 
                               header = TRUE
                               )
                    )

b.pow.imp <- as.numeric(
                        read.delim(
                                   "beta_a10_m1.txt", 
                                   header = TRUE
                                   )
                        )

b.dep <- as.numeric(
                    read.delim(
                               "beta_a9_m5.txt", 
                               header = TRUE
                               )
                    )

b.dep.imp <- as.numeric(
                        read.delim(
                                   "beta_a10_m5.txt", 
                                   header = TRUE
                                   )
                        )

#############
# Figure A5 #
#############

CorrPredPlot(
             sanctions$success,
             x.pow,
             x.dep,
             b.pow,
             b.dep
             )

#############
# Figure A6 #
#############

CorrPredPlot(
             sanctions$success,
             x.pow,
             x.dep,
             b.pow.imp,
             b.dep.imp
             )

##############################
# Compare elasticity and RCA #
##############################

pow.elas <- read.dta13("power_elasticity.dta")
                        # Read in a data set with elasticity and
                        # RCA data combined

pow.el <- data.frame(
                     country   = pow.elas$country,
                     commodity = pow.elas$category,
                     rca       = pow.elas$rca_scaled,
                     elas      = pow.elas$elasticity
                     )

#############
# Figure A7 #
#############

ggplot(
       pow.el,
       aes(
           x = rca,
           y = elas
           )
       ) +
  geom_point(
             size   = 1.25,
             alpha  = 0.5,
             colour = "navyblue"
             ) +
  ylab("Export Elasticity") +
  xlab("Relative Revealed Comparative Advantage") +
  theme_bw()

